
#Replication file for "Ought it Audit," Anson & Kane, JEPS 2022

#1 package requirements

require(tidyverse)
require(car)
require(ggthemes)
require(readxl)
require(sjPlot)
require(stargazer)
library(psych)
require(MASS)
require(gridExtra)
require(pwr)
require(haven)
require(fastDummies)


#2 data read-in

dat <- read.csv("jeps_replication.csv")


#3 Modeling

mod1 <- lm(data=dat, IRS_dv~as.factor(IRSExpConds))
mod2 <- lm(data=dat, IRS_DV2_ReceiveInfo_binary~as.factor(IRSExpConds))
mod3 <- lm(data=dat, IRS_DV3_WriteLetter_binary~as.factor(IRSExpConds))


#4 Replication of Fig. 1

modelLabels <- c("Increase IRS Funding",
                 "Obtain IRS Information",
                 "Representative Contact"
)


modelCoefs <- c(mod1$coefficients[2],
                mod2$coefficients[2],
                mod3$coefficients[2]
)

modelLow95 <- c(confint(mod1,level=0.95)[2,1],
                confint(mod2,level=0.95)[2,1],
                confint(mod3,level=0.95)[2,1]
) 

modelLow90 <- c(confint(mod1,level=0.90)[2,1],
                confint(mod2,level=0.90)[2,1],
                confint(mod3,level=0.90)[2,1]
)

modelHi95 <- c(confint(mod1,level=0.95)[2,2],
               confint(mod2,level=0.95)[2,2],
               confint(mod3,level=0.95)[2,2]
) 

modelHi90 <- c(confint(mod1,level=0.90)[2,2],
               confint(mod2,level=0.90)[2,2],
               confint(mod3,level=0.90)[2,2]
)

grafMat <- data.frame(modelLabels,modelCoefs,modelHi90,
                      modelHi95,modelLow90,modelLow95)

grafMat$xlabs <- modelLabels


graf1 <- grafMat %>%
  ggplot() +
  geom_point(aes(y=reorder(xlabs,desc(xlabs)),x=modelCoefs),size=2, color="Dark Green")+
  geom_errorbarh(aes(y=reorder(xlabs,desc(xlabs)),xmax=modelHi90,xmin=modelLow90),height=0.0,size=1.1, color="Dark Green")+
  geom_errorbarh(aes(y=reorder(xlabs,desc(xlabs)),xmax=modelHi95,xmin=modelLow95),height=0.1, color="Dark Green")+
  geom_vline(aes(xintercept=0),linetype=2)+
  ggtitle("Support for Increased IRS Funding")+
  xlab("Model Coefficient")+
  ylab("")+
  theme_few()+
  theme(panel.background = element_rect(fill='grey95'))+
  theme(text=element_text(size=10))+
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  scale_color_manual(values=c("dark blue","dark red"))+
  theme(legend.position="none")


graf1




#5 Replication of Fig. 2


mod1 <- lm(data=dat[dat$D0R1==1,], IRS_dv~as.factor(IRSExpConds))
mod1b <- lm(data=dat[dat$D0R1==0,], IRS_dv~as.factor(IRSExpConds))
mod2 <- lm(data=dat[dat$D0R1==1,], IRS_DV2_ReceiveInfo_binary~as.factor(IRSExpConds))
mod2b <- lm(data=dat[dat$D0R1==0,], IRS_DV2_ReceiveInfo_binary~as.factor(IRSExpConds))
mod3 <- lm(data=dat[dat$D0R1==1,], IRS_DV3_WriteLetter_binary~as.factor(IRSExpConds))
mod3b <- lm(data=dat[dat$D0R1==0,], IRS_DV3_WriteLetter_binary~as.factor(IRSExpConds))


modelLabels <- c("Increase IRS Funding",
                 "Increase IRS Funding",
                 "Obtain IRS Information",
                 "Obtain IRS Information",
                 "Representative Contact",
                 "Representative Contact"
)



modelCoefs <- c(mod1$coefficients[2],mod1b$coefficients[2],
                mod2$coefficients[2],mod2b$coefficients[2],
                mod3$coefficients[2],mod3b$coefficients[2]
)

modelLow95 <- c(confint(mod1,level=0.95)[2,1],
                confint(mod1b,level=0.95)[2,1],
                confint(mod2,level=0.95)[2,1],
                confint(mod2b,level=0.95)[2,1],
                confint(mod3,level=0.95)[2,1],
                confint(mod3b,level=0.95)[2,1]
) 

modelLow90 <- c(confint(mod1,level=0.90)[2,1],
                confint(mod1b,level=0.90)[2,1],
                confint(mod2,level=0.90)[2,1],
                confint(mod2b,level=0.90)[2,1],
                confint(mod3,level=0.90)[2,1],
                confint(mod3b,level=0.90)[2,1]
)

modelHi95 <- c(confint(mod1,level=0.95)[2,2],
                confint(mod1b,level=0.95)[2,2],
                confint(mod2,level=0.95)[2,2],
                confint(mod2b,level=0.95)[2,2],
                confint(mod3,level=0.95)[2,2],
                confint(mod3b,level=0.95)[2,2]
) 

modelHi90 <- c(confint(mod1,level=0.90)[2,2],
                confint(mod1b,level=0.90)[2,2],
                confint(mod2,level=0.90)[2,2],
                confint(mod2b,level=0.90)[2,2],
                confint(mod3,level=0.90)[2,2],
                confint(mod3b,level=0.90)[2,2]
)

grafMat <- data.frame(modelLabels,modelCoefs,modelHi90,
                      modelHi95,modelLow90,modelLow95)


grafMat$partyLabels <- c("Republicans","Democrats",
                         "Republicans","Democrats",
                         "Republicans","Democrats")


grafMat$row <- c(1,1,2,2,3,3)

grafMat$xlabs <- modelLabels


graf2 <- grafMat %>%
  ggplot() +
  geom_point(aes(y=partyLabels,x=modelCoefs),size=2, color="Dark Green")+
  geom_errorbarh(aes(y=partyLabels,xmax=modelHi90,xmin=modelLow90),height=0.0,size=1.1, color="Dark Green")+
  geom_errorbarh(aes(y=partyLabels,xmax=modelHi95,xmin=modelLow95),height=0.1, color="Dark Green")+
  geom_vline(aes(xintercept=0),linetype=2)+
  ggtitle("Support for Increased IRS Funding by Party ID")+
  xlab("Model Coefficient")+
  ylab("")+
  theme_few()+
  theme(panel.background = element_rect(fill='grey95'))+
  theme(text=element_text(size=10))+
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  scale_color_manual(values=c("dark blue","dark red"))+
  theme(legend.position="none")+
  facet_wrap(~xlabs,ncol=1)


graf2

#6 Replication of Fig. 3
##First row of Fig. 3

dat$IRSExpConds2 <- as.factor(dat$IRSExpConds)
dat$IRSExpConds2 <- relevel(dat$IRSExpConds2, ref="1")

mod1 <- lm(data=dat[dat$D0R1==1,], IRS_dv~as.factor(IRSExpConds2))
mod1b <- lm(data=dat[dat$D0R1==0,], IRS_dv~as.factor(IRSExpConds2))

modelLabels <- c("Deficit Framing",
                 "Deficit Framing",
                 "Inequality Framing",
                 "Inequality Framing"
)


modelCoefs <- c(mod1$coefficients[3],mod1b$coefficients[3],
                mod1$coefficients[4],mod1b$coefficients[4]
)

modelLow95 <- c(confint(mod1,level=0.95)[3,1],
                confint(mod1b,level=0.95)[3,1],
                confint(mod1,level=0.95)[4,1],
                confint(mod1b,level=0.95)[4,1]
) 

modelLow90 <- c(confint(mod1,level=0.90)[3,1],
                confint(mod1b,level=0.90)[3,1],
                confint(mod1,level=0.90)[4,1],
                confint(mod1b,level=0.90)[4,1]
) 

modelHi95 <- c(confint(mod1,level=0.95)[3,2],
                confint(mod1b,level=0.95)[3,2],
                confint(mod1,level=0.95)[4,2],
                confint(mod1b,level=0.95)[4,2]
) 

modelHi90 <- c(confint(mod1,level=0.90)[3,2],
                confint(mod1b,level=0.90)[3,2],
                confint(mod1,level=0.90)[4,2],
                confint(mod1b,level=0.90)[4,2]
) 

grafMat <- data.frame(modelLabels,modelCoefs,modelHi90,
                      modelHi95,modelLow90,modelLow95)


grafMat$partyLabels <- c("Republicans","Democrats",
                         "Republicans","Democrats")


grafMat$row <- c(1,1,2,2)

grafMat$xlabs <- modelLabels


graf3 <- grafMat %>%
  ggplot() +
  geom_point(aes(y=reorder(xlabs,desc(xlabs)),x=modelCoefs),size=2, color="Dark Green")+
  geom_errorbarh(aes(y=reorder(xlabs,desc(xlabs)),xmax=modelHi90,xmin=modelLow90),height=0.0,size=1.1, color="Dark Green")+
  geom_errorbarh(aes(y=reorder(xlabs,desc(xlabs)),xmax=modelHi95,xmin=modelLow95),height=0.1, color="Dark Green")+
  geom_vline(aes(xintercept=0),linetype=2)+
  ggtitle("Support for Increased IRS Funding by Party ID")+
  xlab("")+
  ylab("")+
  scale_x_continuous(limits=c(-0.25,0.3),breaks=c(-0.2,-0.1,0,0.1,0.2,0.3))+
  theme_few()+
  theme(panel.background = element_rect(fill='grey95'))+
  theme(text=element_text(size=10))+
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  scale_color_manual(values=c("dark blue","dark red"))+
  theme(legend.position="none")+
  facet_wrap(~reorder(partyLabels,desc(partyLabels)),ncol=2)



##Second row of Fig. 3


mod1 <- lm(data=dat[dat$D0R1==1,], IRS_DV2_ReceiveInfo_binary~IRSExpConds2)
mod1b <- lm(data=dat[dat$D0R1==0,], IRS_DV2_ReceiveInfo_binary~IRSExpConds2)

modelLabels <- c("Deficit Framing",
                 "Deficit Framing",
                 "Inequality Framing",
                 "Inequality Framing"
)

modelCoefs <- c(mod1$coefficients[3],mod1b$coefficients[3],
                mod1$coefficients[4],mod1b$coefficients[4]
)

modelLow95 <- c(confint(mod1,level=0.95)[3,1],
                confint(mod1b,level=0.95)[3,1],
                confint(mod1,level=0.95)[4,1],
                confint(mod1b,level=0.95)[4,1]
) 

modelLow90 <- c(confint(mod1,level=0.90)[3,1],
                confint(mod1b,level=0.90)[3,1],
                confint(mod1,level=0.90)[4,1],
                confint(mod1b,level=0.90)[4,1]
) 

modelHi95 <- c(confint(mod1,level=0.95)[3,2],
               confint(mod1b,level=0.95)[3,2],
               confint(mod1,level=0.95)[4,2],
               confint(mod1b,level=0.95)[4,2]
) 

modelHi90 <- c(confint(mod1,level=0.90)[3,2],
               confint(mod1b,level=0.90)[3,2],
               confint(mod1,level=0.90)[4,2],
               confint(mod1b,level=0.90)[4,2]
) 

grafMat <- data.frame(modelLabels,modelCoefs,modelHi90,
                      modelHi95,modelLow90,modelLow95)


grafMat$partyLabels <- c("Republicans","Democrats",
                         "Republicans","Democrats")


grafMat$row <- c(1,1,2,2)

grafMat$xlabs <- modelLabels


graf4 <- grafMat %>%
  ggplot() +
  geom_point(aes(y=reorder(xlabs,desc(xlabs)),x=modelCoefs),size=2, color="Dark Green")+
  geom_errorbarh(aes(y=reorder(xlabs,desc(xlabs)),xmax=modelHi90,xmin=modelLow90),height=0.0,size=1.1, color="Dark Green")+
  geom_errorbarh(aes(y=reorder(xlabs,desc(xlabs)),xmax=modelHi95,xmin=modelLow95),height=0.1, color="Dark Green")+
  geom_vline(aes(xintercept=0),linetype=2)+
  ggtitle("Willingness to Receive IRS Info by Party ID")+
  xlab("")+
  ylab("")+
  scale_x_continuous(limits=c(-0.25,0.3),breaks=c(-0.2,-0.1,0,0.1,0.2,0.3))+
  theme_few()+
  theme(panel.background = element_rect(fill='grey95'))+
  theme(text=element_text(size=10))+
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  scale_color_manual(values=c("dark blue","dark red"))+
  theme(legend.position="none")+
  facet_wrap(~reorder(partyLabels,desc(partyLabels)),ncol=2)


##Bottom row of Fig. 3


mod1 <- lm(data=dat[dat$D0R1==1,], IRS_DV3_WriteLetter_binary~as.factor(IRSExpConds2))
mod1b <- lm(data=dat[dat$D0R1==0,], IRS_DV3_WriteLetter_binary~as.factor(IRSExpConds2))

modelLabels <- c("Deficit Framing",
                 "Deficit Framing",
                 "Inequality Framing",
                 "Inequality Framing"
)


modelCoefs <- c(mod1$coefficients[3],mod1b$coefficients[3],
                mod1$coefficients[4],mod1b$coefficients[4]
)

modelLow95 <- c(confint(mod1,level=0.95)[3,1],
                confint(mod1b,level=0.95)[3,1],
                confint(mod1,level=0.95)[4,1],
                confint(mod1b,level=0.95)[4,1]
) 

modelLow90 <- c(confint(mod1,level=0.90)[3,1],
                confint(mod1b,level=0.90)[3,1],
                confint(mod1,level=0.90)[4,1],
                confint(mod1b,level=0.90)[4,1]
) 

modelHi95 <- c(confint(mod1,level=0.95)[3,2],
               confint(mod1b,level=0.95)[3,2],
               confint(mod1,level=0.95)[4,2],
               confint(mod1b,level=0.95)[4,2]
) 

modelHi90 <- c(confint(mod1,level=0.90)[3,2],
               confint(mod1b,level=0.90)[3,2],
               confint(mod1,level=0.90)[4,2],
               confint(mod1b,level=0.90)[4,2]
) 

grafMat <- data.frame(modelLabels,modelCoefs,modelHi90,
                      modelHi95,modelLow90,modelLow95)


grafMat$partyLabels <- c("Republicans","Democrats",
                         "Republicans","Democrats")


grafMat$row <- c(1,1,2,2)

grafMat$xlabs <- modelLabels


graf5 <- grafMat %>%
  ggplot() +
  geom_point(aes(y=reorder(xlabs,desc(xlabs)),x=modelCoefs),size=2, color="Dark Green")+
  geom_errorbarh(aes(y=reorder(xlabs,desc(xlabs)),xmax=modelHi90,xmin=modelLow90),height=0.0,size=1.1, color="Dark Green")+
  geom_errorbarh(aes(y=reorder(xlabs,desc(xlabs)),xmax=modelHi95,xmin=modelLow95),height=0.1, color="Dark Green")+
  geom_vline(aes(xintercept=0),linetype=2)+
  ggtitle("Willingness to Contact Representatives by Party ID")+
  xlab("Model Coefficient")+
  ylab("")+
  scale_x_continuous(limits=c(-0.25,0.3),breaks=c(-0.2,-0.1,0,0.1,0.2,0.3))+
  theme_few()+
  theme(panel.background = element_rect(fill='grey95'))+
  theme(text=element_text(size=10))+
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  scale_color_manual(values=c("dark blue","dark red"))+
  theme(legend.position="none")+
  facet_wrap(~reorder(partyLabels,desc(partyLabels)),ncol=2)



##Assembly of Fig. 3

grid.arrange(graf3,graf4,graf5,ncol=1)





#SUPPLEMENTARY INFORMATION

dat$LibCon <- car::recode(dat$LibCon,"-99=NA")
dat$ideo <- car::recode(dat$LibCon, "1:2=1;3=2;4:5=3")
#S1 Replication of Table S1

ethnicity <- dummy_cols(dat$ethnicity)
ethnicity <- ethnicity[,2:3]
colnames(ethnicity) <- c("white","black")
pid <- dummy_cols(dat$D0I1R2)
pid <- pid[,2:4]
colnames(pid) <- c("Democrat","Independent","Republican")
ideo <- dummy_cols(dat$ideo)
ideo <- ideo[,2:4]
colnames(ideo) <- c("Liberal","Moderate","Conservative")
region <- dummy_cols(dat$region)
region <- region[,2:5]
colnames(region) <- c("Northeast","South","Midwest","West")

dat2 <- data.frame(dat$hhi,dat$age,dat$gender,ethnicity,dat$hispanic,
                   pid,ideo,region)

dat2$dat.gender <- dat2$dat.gender-1
dat2$dat.hispanic <- car::recode(dat2$dat.hispanic,"1=0;else=1")

des <- psych::describe(dat2)

#based on this info., and U.S. Census & 2021 Gallup data, construct a
#table containing each median value or proportion:

describe <- data.frame(colnames(dat2),des$mean)
describe$des.mean <- describe$des.mean*100
describe$des.mean[1:2] <- c("$35k-39k",des$median[2])

describe$benchmarks <- c("$67k","38.1","51%","76.3%","13.4%","18.5%",
                         "42%","11%","47%","","","","20%","34%","20%","26%")

write.table(describe,"table_s1.csv",sep=",",col.names = T, row.names = T)


#SI Section 2 


#Table S2. manipulation checks

mod1 <- lm(data=dat,IRS_fmc~as.factor(IRSExpConds))
mod2 <- lm(data=dat,IRS_smc~as.factor(IRSExpConds))
mod3 <- lm(data=dat,IRS_smc2~as.factor(IRSExpConds))
mod4 <- lm(data=dat,IRS_smc3~as.factor(IRSExpConds))

stargazer(mod1,mod2,mod3,mod4,
          type="html",out="Table_S2.html",
          star.cutoffs = c(0.1,0.05,0.01))


#Table S3. Manipulation Check by Party

mod1 <- lm(data=dat[dat$D0R1==0,],IRS_fmc~as.factor(IRSExpConds))
mod2 <- lm(data=dat[dat$D0R1==0,],IRS_smc2~as.factor(IRSExpConds))
mod3 <- lm(data=dat[dat$D0R1==0,],IRS_smc3~as.factor(IRSExpConds))
mod1b <- lm(data=dat[dat$D0R1==1,],IRS_fmc~as.factor(IRSExpConds))
mod2b <- lm(data=dat[dat$D0R1==1,],IRS_smc2~as.factor(IRSExpConds))
mod3b <- lm(data=dat[dat$D0R1==1,],IRS_smc3~as.factor(IRSExpConds))

stargazer(mod1,mod2,mod3,mod1b,mod2b,mod3b,
         type="html",out="Table_S3.html",
         star.cutoffs = c(0.1,0.05,0.01))




######SI Section 3

#Replication of Table S4. Tabular OLS Regression Corresponding to Fig. 1

mod1 <- lm(data=dat, IRS_dv~as.factor(IRSExpConds))
mod2 <- lm(data=dat, IRS_DV2_ReceiveInfo_binary~as.factor(IRSExpConds))
mod3 <- lm(data=dat, IRS_DV3_WriteLetter_binary~as.factor(IRSExpConds))

stargazer(mod1,mod2,mod3,type="html",out="Table_S4.html",
         star.cutoffs = c(0.1,0.05,0.01))




#Replication of Table S5. Tabular OLS Regression Corresponding to Fig. 2 and Fig. 3

dat$IRSExpConds2 <- as.factor(dat$IRSExpConds)
dat$IRSExpConds2 <- relevel(dat$IRSExpConds2, ref = "1")

mod1 <- lm(data=dat[dat$D0R1==1,], IRS_dv~as.factor(IRSExpConds2))
mod1b <- lm(data=dat[dat$D0R1==0,], IRS_dv~as.factor(IRSExpConds2))
mod2 <- lm(data=dat[dat$D0R1==1,], IRS_DV2_ReceiveInfo_binary~as.factor(IRSExpConds2))
mod2b <- lm(data=dat[dat$D0R1==0,], IRS_DV2_ReceiveInfo_binary~as.factor(IRSExpConds2))
mod3 <- lm(data=dat[dat$D0R1==1,], IRS_DV3_WriteLetter_binary~as.factor(IRSExpConds2))
mod3b <- lm(data=dat[dat$D0R1==0,], IRS_DV3_WriteLetter_binary~as.factor(IRSExpConds2))

stargazer(mod1,mod1b,mod2,mod2b,mod3,mod3b,type="html",out="Table_S5.html",
         star.cutoffs = c(0.1,0.05,0.01))

#Replication of Table S6. Alternative Model Specifications


#Need to reverse the coding of the three-item DV2 and DV3 scales; not yet
#cleaned in the data
dat$IRS_dv2 <- car::recode(dat$IRS_dv2, "1=3;2=2;3=1")
dat$IRS_dv3 <- car::recode(dat$IRS_dv3, "1=3;2=2;3=1")

modO1 <- lm(data=dat, IRS_DV2_ReceiveInfo_binary~as.factor(IRSExpConds))
modO2 <- lm(data=dat, IRS_DV3_WriteLetter_binary~as.factor(IRSExpConds))
modB1 <- glm(data=dat, IRS_DV2_ReceiveInfo_binary~as.factor(IRSExpConds), family=binomial(link="logit"))
modB2 <- glm(data=dat, IRS_DV3_WriteLetter_binary~as.factor(IRSExpConds), family=binomial(link="logit"))
modX1 <- lm(data=dat, IRS_dv2~as.factor(IRSExpConds))
modX2 <- lm(data=dat, IRS_dv3~as.factor(IRSExpConds))
modZ1 <- polr(data=dat, as.factor(IRS_dv2)~as.factor(IRSExpConds))
modZ2 <- polr(data=dat, as.factor(IRS_dv3)~as.factor(IRSExpConds))

stargazer(modO1,modO2,modB1,modB2,modX1,modX2,
         modZ1,modZ2,type="html",out="Table_S6.html",
         star.cutoffs = c(0.1,0.05,0.01))




#Table S7. Treatment Effects and PID Strength

#Republicans
modP1 <- lm(data=dat, IRS_dv ~ as.factor(IRSExpConds2) + RepIdentityclean +
              as.factor(IRSExpConds2)*RepIdentityclean)
#Democrats
modP2 <- lm(data=dat, IRS_dv ~ as.factor(IRSExpConds2) + DemIdentityclean +
              as.factor(IRSExpConds2)*DemIdentityclean)

stargazer(modP1,modP2,
         type="html",out="Table_S7.html",
         star.cutoffs = c(0.1,0.05,0.01))




#############SI Section 5: MTurk Experimental Pre-Test

#read in data
dat <- read.csv("jeps_replication_mturk_pretest.csv")



##Fig. S1: Results from experimental pre-test

mod1 <- lm(data=dat, IRS_dv~as.factor(IRSconditions))
mod2 <- lm(data=dat, IRS_dv2_dichot~as.factor(IRSconditions))
mod3 <- lm(data=dat, IRS_dv2_dichot~as.factor(IRSconditions))


modelLabels <- c("Increase IRS Funding",
                 "Obtain IRS Information",
                 "Representative Contact"
)


modelCoefs <- c(mod1$coefficients[2],
                mod2$coefficients[2],
                mod3$coefficients[2]
)

modelLow95 <- c(confint(mod1,level=0.95)[2,1],
                confint(mod2,level=0.95)[2,1],
                confint(mod3,level=0.95)[2,1]
) 

modelLow90 <- c(confint(mod1,level=0.90)[2,1],
                confint(mod2,level=0.90)[2,1],
                confint(mod3,level=0.90)[2,1]
)

modelHi95 <- c(confint(mod1,level=0.95)[2,2],
               confint(mod2,level=0.95)[2,2],
               confint(mod3,level=0.95)[2,2]
) 

modelHi90 <- c(confint(mod1,level=0.90)[2,2],
               confint(mod2,level=0.90)[2,2],
               confint(mod3,level=0.90)[2,2]
)

grafMat <- data.frame(modelLabels,modelCoefs,modelHi90,
                      modelHi95,modelLow90,modelLow95)

grafMat$xlabs <- modelLabels

graf1 <- grafMat %>%
  ggplot() +
  geom_point(aes(y=reorder(xlabs,desc(xlabs)),x=modelCoefs),size=2, color="Dark Green")+
  geom_errorbarh(aes(y=reorder(xlabs,desc(xlabs)),xmax=modelHi90,xmin=modelLow90),height=0.0,size=1.1, color="Dark Green")+
  geom_errorbarh(aes(y=reorder(xlabs,desc(xlabs)),xmax=modelHi95,xmin=modelLow95),height=0.1, color="Dark Green")+
  geom_vline(aes(xintercept=0),linetype=2)+
  ggtitle("T1 vs. Control")+
  xlab("Model Coefficient")+
  ylab("")+
  theme_few()+
  theme(panel.background = element_rect(fill='grey95'))+
  theme(text=element_text(size=10))+
  theme(plot.title = element_text(hjust = 0.5,size=10))+
  scale_color_manual(values=c("dark blue","dark red"))+
  theme(legend.position="none")

#create graphic
graf1



##Table S9:  Interactions Between Party and
#Deficit- vs. Inequality-Reduction Frames, Amazon MTurk Sample

mod1 <- lm(data=dat[dat$IRSconditions==2|dat$IRSconditions==3,], IRS_dv_clean~as.factor(IRSconditions)+D0R1 + as.factor(IRSconditions)*D0R1)
mod2 <- polr(data=dat[dat$IRSconditions==2|dat$IRSconditions==3,], as.factor(IRS_dv2)~as.factor(IRSconditions)+D0R1 + as.factor(IRSconditions)*D0R1)
mod3 <- polr(data=dat[dat$IRSconditions==2|dat$IRSconditions==3,], as.factor(IRS_dv3)~as.factor(IRSconditions)+D0R1 + as.factor(IRSconditions)*D0R1)

stargazer(mod1,mod2,mod3,
          type="html",out="Table_S9.html",
          star.cutoffs = c(0.1,0.05,0.01))





######SI section 6: Power Analysis

#We averaged the R^2 values produced in the MTurk study's comparison of control
#vs. T1, control vs. T2, and control vs. T3 in the full sample: 

(0.0338 + 0.0263 + 0.0282)/3

#The result was an average R^2 of 0.0294

#We estimate an F2 value of 0.03:
f2val <- (0.0294/(1-0.0294))

#We compare Power F2 tests at power levels 0.7, 0.8, and 0.9
pwr.f2.test(u = 5, f2 = f2val, sig.level = 0.05, power = 0.7)
pwr.f2.test(u = 5, f2 = f2val, sig.level = 0.05, power = 0.8)
pwr.f2.test(u = 5, f2 = f2val, sig.level = 0.05, power = 0.9)

#We also performed a power F2 test with the specified N = 1,050.
pwr.f2.test(u = 5, v = 1050, f2 = f2val, sig.level = 0.05)





###########SI Section 7: Results from historical surveys

#Determining the partisan perceptions of the IRS, Pew 2015 study

dat <- read_spss("pew2015.por")

#This is Pew's 2015 September Political Survey

table(dat$PARTY)
dat$PARTY <- as.numeric(dat$PARTY)
dat$pid3 <- car::recode(dat$PARTY,"1='Republican';2='Democrat';else=NA")
dat$PARTYLN <- as.numeric(dat$PARTYLN)

dat$pid3[dat$PARTYLN==1] <- "Republican"
dat$pid3[dat$PARTYLN==2] <- "Democrat"

dat$irs <- as.numeric(dat$Q13OF2)
dat$irs <- car::recode(dat$irs, "1=4;2=3;3=2;4=1;else=NA")


#Summarize attitudes by party
dat2 <- dat %>%
  dplyr::select(irs,pid3) %>%
  na.omit() %>%
  group_by(pid3) %>%
  summarize(irsMean = mean(irs),
            irsSD = sd(irs),
            n = n())

dat2$year <- "2015"


#Repeat the process with Pew data from 2018
dat <- read.csv("pew2018.csv")

#This is the Pew Supreme Court and Tariffs July 2018 Survey

table(dat$party)
dat$party <- as.numeric(dat$party)
dat$pid3 <- car::recode(dat$party,"2='Republican';1='Democrat';else=NA")
dat$partyln <- as.numeric(dat$partyln)

dat$pid3[dat$PARTYLN==4] <- "Republican"
dat$pid3[dat$PARTYLN==2] <- "Democrat"

dat$irs <- as.numeric(dat$forq1d)
dat$irs <- car::recode(dat$irs, "7=1;3=2;2=3;6=4;else=NA")


dat3 <- dat %>%
  dplyr::select(irs,pid3) %>%
  na.omit() %>%
  group_by(pid3) %>%
  summarize(irsMean = mean(irs),
            irsSD = sd(irs),
            n = n())

dat3$year <- "2018"



#Bind the two surveys together to perform a single graphical analysis

dat4 <- rbind(dat2,dat3)


#Calculate 95% confidence intervals for the figure

dat4$hi <- dat4$irsMean + 1.96*dat4$irsSD/sqrt(dat4$n)
dat4$lo <- dat4$irsMean - 1.96*dat4$irsSD/sqrt(dat4$n)

#Create the graphic
dat4 %>%
  ggplot() +
  geom_bar(aes(x=pid3,y=irsMean,fill=pid3),stat="identity")+
  geom_errorbar(aes(x=pid3,ymin=lo,ymax=hi),width=0.1)+
  facet_wrap(~year,ncol=2)+
  ylab("Average IRS Support (1:4)")+
  xlab("Party ID (Includes Leaners)")+
  scale_fill_grey()+
  theme_few()+
  theme(legend.position = "none")




